clear all
close all
path1=['C:\Users\Tom�s\Desktop\temp\'];
% range=10;  %std times range for plotting
% initialFile=17
% trials=10
fileName=' '

dir_to_search=[path1]
formatpattern = fullfile(dir_to_search, '/*.ibw');
dInfoNS = dir(formatpattern);
[~, idx] =sort([dInfoNS.datenum]); % sort structure indices by date of creation, to avoid 10 coming before 2.

dinfo=dInfoNS(idx); % turn the not sorted structure to a sorted one.

trials=numel(dinfo)
%% this section just averages Igor waves contained in a directory, including subdirectories?

for  i=1:numel(dinfo)
   
    fileName = [dinfo(i).folder '\' dinfo(i).name];
    Ibw= IBWread(fileName); %49
    dataMat(:,i)=Ibw.y;
    
end
dataMat=medfilt1(dataMat,5);

meanTrace=((mean(dataMat,2)));
t=(1:numel(meanTrace))*Ibw.dx*1000;
%%
%find stim artifact peak location
 [pks locs]=findpeaks(meanTrace,'minPeakHeight',200,'minpeakDistance',200);
locs=locs-5
% locs=[100000 102000 104000 106000 108000]'
locs=[2000 2400 2800];
%remove stim artifact from raw traces
for i=1:numel(dinfo)
    for j=1:numel(locs)
    dataMat([locs(j)-20:locs(j)+20],i)=nan;
%    dataMat(dataMat>15)=NaN;
    dataMat(:,i)=dataMat(:,i)-nanmean(dataMat([1:2e3],i));
    end
end

dataMat=dataMat*1e12;
dataMat(dataMat>10)=nan;
meanTrace=((nanmean(dataMat,2))); % update average after stim artifact removal

meanTrace(meanTrace>8)=nan;
%% Analyze lag for first epsc
pkprom=20 ; %conservative threshold for prominence

%use locs from before, that represent indices of stim artifact peak.
epscLag=[];
epscSize=[];
for i=1:numel(dinfo)
    for j=1:numel(locs); % should be 5 for these trains
    [pks2 locs2]=findpeaks(dataMat([locs(j):end],i)*-1,'minPeakprominence',pkprom,'minpeakDistance',10);
    if ~isempty(locs2)
    epscLag=[epscLag;min(locs2)]; %append the first epsc peak location relative to onset of stimulation.
    epscSize=[epscSize;pks2(locs2==min(locs2))]; %take the amplitude of that epsc.
    else
        epscLag=[epscLag;NaN];
        epscSize=[epscSize;NaN];
    end 
    end
end

epscLag(epscLag>(20/(Ibw.dx*1000)))=NaN;
epscSize(epscLag>(20/(Ibw.dx*1000)))=NaN;  %delete trials where there is no epsc.
d{1}=epscLag*Ibw.dx*1000;
% %% Make a pretty figure :)
% figure;
% raincloud_plot('X', d{1}, 'box_on', 1, 'box_dodge', 1, 'box_dodge_amount',...
% 0, 'dot_dodge_amount', .3, 'color', rgb('deepskyblue'), 'cloud_edge_col', rgb('deepskyblue'));
% figure
% raincloud_plot('X', d{1}, 'box_on', 1);


figure;
hold on
for i=1:numel(dinfo)
    plot(t,dataMat(:,i),'color',[0 0 0 2/numel(dinfo)],'linewidth',2);
end
plot(t,meanTrace,'color',(rgb('red')),'linewidth',2)
xlim([990 1100])
box off
axis off
set(gcf,'color','white');

for i=1:numel(dinfo)
plot((locs+epscLag((i-1)*numel(locs)+1:(i-1)*numel(locs)+numel(locs)))*Ibw.dx*1000,20,'ko')
end

 sb=scalebar('XUnit','ms','YUnit','pA','YLen',50)
title(Ibw.WaveNotes(6:16))  

%%
epscLag(epscLag>(20/(Ibw.dx*1000)))=NaN; %delete epscLag values where failures ocurr.

% %
% figure;
% plot(epscLag*Ibw.dx*1000,1,'ko','markersize',10);
% hold on
% vline(nanmean(epscLag*Ibw.dx*1e3),'r--');
% vline(nanmedian(epscLag*Ibw.dx*1e3),'b--');
% set(gcf,'color','white');
% set(gca,'FontSize', 18)
% box off
% xlim([0 10])
% yticklabels({})
% xlabel('epsc lag (ms)')
% 

clipboard('copy', epscLag)
fclose all


%% d
figure
h=histogram(d{1},[0:0.2:10],'facecolor',rgb('deepskyblue'))
xlim([0 10])
y=ylim;
hold on
plot(d{1}, -y(2)/3+(rand(numel(d{1}),1)*(y(2)/6)),'o','color',[0.4 0.4 0.4 0.5])
box off
[val loc]=max(h.Values);
histpeak=(loc*h.BinWidth)+h.BinLimits(1)-h.BinWidth/2;
vline(histpeak,'--r');
title(Ibw.WaveNotes)





%% section to append new data to mat file with all cells.

% path='C:\Users\tomas\Documents\1. Harvard\2017 - Spring\Regehr lab rotation\Data summary'
% load([path  '\MossyFiberStimV2.mat'])
% numb=nansum(MLILag,2);
% nMLI=sum(numb>0)
% 
% newData=epscLag*Ibw.dx*1000;
% newData2=epscSize;
% 
% MLILag(nMLI+1,1:numel(newData))=newData;
% MLIAmplitude(nMLI+1,1:numel(newData))=newData2;
% MLInTrials(nMLI+1)=trials*numel(locs);
% 
% 
% save([path  '\MossyFiberStimV2.mat'], 'CCsLag','MLILag','CCsAmplitude','MLIAmplitude','CCsnTrials','MLInTrials')

